########
# Replicates conjoint analysis results
#######

#
rm(list=ls())


# load libraries
packs = c('lmtest', 'tidyverse', 'hrbrthemes','interplot',
          'broom', "here", "sparsereg",
          'stargazer')

# load custom functions
source("repFile/r/LoadPkg.R")
source("repFile/r/vcovCluster.R")

#
loadPkg(packs)

# read data
load(here("repFile", "data" , "conjoint-nodup.rda"))
df = conjoint_nodup
rm(conjoint_nodup)



# ivs
covars_fac = c('tierras_fac', 'elections_fac', 'reint_fac', 'just_fac', 
               'repar_fac', 'drugs_fac')




# Conjoint model using factor attributes ----------------------------------


m1 = lm(chosen_dummy ~ just_fac +
          tierras_fac + elections_fac +
          reint_fac + repar_fac +
          drugs_fac, data = df)

summary(m1)


var_replace = c('Intercept', 
                'jail for HR violators', 
                'jail for all FARC', 
                'small transfer', 'large transfer', 
                'FARC competes', 'FARC competes + seats',
                '90\\% min. wage for demob.', 
                '200\\% min. wage for demob.',
                'FARC ask forgiveness', 'FARC pay victims',
                'Manual erradication', 'Substitution program')


# correct standard errors
m1_std = coeftest(m1, vcov. = vcovCluster(model = m1, 
                                          cluster = df$case_id))


# standard tables
stargazer(m1_std, 
          title = 'Average marginal component effects. Standard errors clustered by respondent.', 
          label = 'stdreg',
          style = 'ajps',
          covariate.labels = var_replace,
          intercept.bottom = F, 
          out = here("repFile", "paper", "figures", "amces.tex"))


# standard tables
stargazer(m1_std, 
          title = 'Average marginal component effects. Standard errors clustered by respondent.', 
          label = 'stdreg',
          style = 'ajps',
          covariate.labels = var_replace,
          intercept.bottom = F, 
          out = here("repFile", "JOP submission", "amces.tex"))



# ols interaction model
m1 = lm(chosen_dummy ~ just_fac*rearing_agg +
          tierras_fac*rearing_agg + elections_fac*rearing_agg +
          reint_fac*rearing_agg + repar_fac*rearing_agg +
          drugs_fac*rearing_agg, data = df)

summary(m1)


var_replace = c('Intercept', 
                'jail for HR violators', 
                'jail for all FARC', 
                'Child-Rearing Score', 
                'small transfer', 'large transfer', 
                'FARC competes', 'FARC competes + seats',
                '90\\% min. wage for demob.', 
                '200\\% min. wage for demob.',
                'FARC ask forgiveness', 'FARC pay victims',
                'Manual erradication', 'Substitution program', 
                'Child-Rearing Score X Jail for HR', 
                'Child-Rearing Score X Jail for all', 
                'Child-Rearing Score X Low transfer', 
                'Child-Rearing Score X High transfer',
                'Child-Rearing Score X FARC competes', 
                'Child-Rearing Score X FARC seats', 
                'Child-Rearing Score X Low wage', 
                'Child-Rearing Score X High wage', 
                'Child-Rearing Score X Apology', 
                'Child-Rearing Score X Material Repar.', 
                'Child-Rearing Score X Manual erad.', 
                'Child-Rearing Score X Subs. program'
)


# correct standard errors
m1_std = coeftest(m1, vcov. = vcovCluster(model = m1, 
                                          cluster = df$case_id))

# stargazer for tables
stargazer(m1_std,
         title = 'Child-Rearing Interaction Models',
         out = here("repFile", "paper", "figures", "auth-models.tex"),
         label = 'auth-models',
         font.size = "scriptsize",
         style = 'io',
         keep.stat = 'n',
         covariate.labels = var_replace, 
         intercept.bottom = F)

stargazer(m1_std,
          title = 'Child-Rearing Interaction Models',
          out = here("repFile", "JOP submission", "auth-models.tex"),
          label = 'auth-models',
          font.size = "scriptsize",
          style = 'io',
          keep.stat = 'n',
          covariate.labels = var_replace, 
          intercept.bottom = F)



# marginal effect plots
mfxJust = interplot(m1, var1 = 'just_fac', 
                    var2 = 'rearing_agg', point = T)$data
mfxRepar = interplot(m1, var1 = 'reint_fac', 
                     var2 = 'rearing_agg', point = T)$data
mfxElec = interplot(m1, var1 = 'elections_fac', 
                    var2 = 'rearing_agg', point = T)$data



# combine for plotting
pDat = 
  rbind(mfxJust, mfxRepar, mfxElec) %>% 
  filter(str_detect(value, '-high'))

ggplot(pDat, aes(x = fake, y = coef1, color = value, shape = value)) + 
  geom_point(size = 4) + 
  geom_linerange(aes(ymin = ub, ymax = lb)) + 
  facet_grid(~value) + 
  theme_ipsum_ps(grid = 'Y', base_size = 12) + 
  scale_color_manual(labels = c('All FARC face jail time', 
                                'Double minimum salary',
                                'Electoral participation + seats'), 
                       values = wesanderson::wes_palette(name = 'Darjeeling1', 
                                                         n = 3)) + 
  scale_shape_manual(labels = c('All FARC face jail time', 
                                'Double minimum salary',
                                'Electoral participation + seats'), 
                     values = c(17, 16, 15)) + 
  theme(legend.position = 'top',
        legend.title = element_blank(), 
        strip.text = element_blank()) + 
  labs(NULL, 
       x = 'Child-rearing score',
       y = 'Effect of provision on Pr(choose agreement)') + 
  scale_y_percent(limits = c(-.20, .20))
ggsave(here("repFile", "paper", "figures", "mfx-auth-inter.pdf"), 
       device = cairo_pdf)
ggsave(here("repFile", "JOP submission", "mfx-auth-inter.eps"), 
       device = cairo_ps)
ggsave(here("repFile", "JOP submission", "mfx-auth-inter.pdf"), 
       device = cairo_pdf)



# Conjoint model (numeric attributes) -------------------------------------

m1 = lm(chosen_dummy ~ just_num*rearing_agg +
          tierras_num*rearing_agg + elections_num*rearing_agg +
          reint_num*rearing_agg + repar_num*rearing_agg +
          drugs_num*rearing_agg, data = df)

summary(m1)


var_replace = c('Intercept', 
                'Retribution', 
                'Child-rearing', 
                'Land-reform', 
                'Elections',
                'Reintegration', 
                'Reparations', 
                'Drug-reform', 
                'Retribution X child-rearing',
                'Land-reform X child-rearing',
                'Elections X child-rearing',
                'Reintegration X child-rearing',
                'Reparations X child-rearing',
                'Drug-reform X child-rearing'
)




# correct standard errors
m1_std = coeftest(m1, vcov. = vcovCluster(model = m1, 
                                          cluster = df$case_id))

# stargazer for tables
stargazer(m1_std,
         title = 'Child-rearing interaction models. Attributes treated as ordinal. Standard errors clustered at participant-level.',
         out = here("repFile", "paper", "figures", "auth-models-num.tex"),
         label = 'auth-models-num',
         style = 'io',
         covariate.labels = var_replace, 
         intercept.bottom = FALSE, 
         add.lines = list(c("N:", nrow(m1$model))))

stargazer(m1_std,
          title = 'Child-rearing interaction models. Attributes treated as ordinal. Standard errors clustered at participant-level.',
          out = here("repFile", "JOP submission", "auth-models-num.tex"),
          label = 'auth-models-num',
          style = 'io',
          covariate.labels = var_replace, 
          intercept.bottom = FALSE, 
          add.lines = list(c("N:", nrow(m1$model))))



# format data for plotting
pDat = 
  tidy(m1) %>% 
  mutate(label = var_replace) %>% 
  mutate(type = ifelse(str_detect(label, 'X'), 
                       'Interaction terms', 'Component terms')) %>% 
  mutate(lo95 = estimate - `std.error`*1.96, 
         hi95 = estimate + `std.error`*1.96) %>% 
  filter(label != "Intercept") %>% 
  #filter(label != "Child-rearing") %>% 
  #filter(type == 'Interaction') %>% 
  mutate(attribute = word(label, 1)) %>% 
  mutate(sig = ifelse(`p.value` < .05, 1, 0)) %>% 
  mutate(type = factor(type, levels = c('Component terms', 
                                        'Interaction terms')))


ggplot(data = pDat, aes(x = attribute, y = estimate, 
                        ymin = lo95, ymax = hi95, 
                        color = type, shape = factor(sig))) + 
  geom_pointrange(size = .8) + 
  coord_flip() + 
  theme_light(base_family = "Roboto Condensed", base_size = 14) + 
  geom_hline(yintercept = 0, lty = 2) + 
  facet_wrap(~type, scales = 'free') + 
  theme(legend.position = 'none', strip.text = element_text(face = "bold"),
        strip.background = element_rect(colour = "black", fill = "black"), 
        panel.grid.major.x = element_blank(), 
        panel.grid.minor.x = element_blank()) + 
  scale_color_manual(values = wesanderson::wes_palette(name = 'Darjeeling1')) + 
  labs(y = 'Coefficient estimates', 
       x = '')
ggsave(here("repFile", "paper", "figures", "auth-inter-coefs.pdf"),
       device = cairo_pdf, width = 10, height = 6)
ggsave(here("repFile", "JOP submission", "auth-inter-coefs.eps"),
       device = cairo_ps, width = 10, height = 6)
ggsave(here("repFile", "JOP submission", "auth-inter-coefs.pdf"),
       device = cairo_pdf, width = 10, height = 6)



# maginal effects plot
mfxJust = 
  interplot(m1, var1 = 'just_num', var2 = 'rearing_agg', point = T)$data %>% 
  mutate(type = 'Retributive Justice')
mfxRepar = 
  interplot(m1, var1 = 'reint_num', var2 = 'rearing_agg', point = T)$data %>% 
  mutate(type = 'Reparations')
mfxElec = 
  interplot(m1, var1 = 'elections_num', var2 = 'rearing_agg', point = T)$data %>% 
  mutate(type = 'Elections')


# combine for plotting
pDat = 
  rbind(mfxJust, mfxRepar, mfxElec)

ggplot(pDat, aes(x = fake, y = coef1, color = type, shape = type)) + 
  geom_point(size = 4) + 
  geom_linerange(aes(ymin = ub, ymax = lb)) + 
  facet_grid(~type) + 
  theme_ipsum_rc(grid = 'Y') + 
  scale_color_manual(values = wesanderson::wes_palette(name = 'Darjeeling1', 
                                                       n = 3)) + 
  theme(legend.position = 'top',
        legend.title = element_blank(), 
        strip.text = element_blank()) + 
  labs(NULL, 
       x = 'Child-rearing score',
       y = 'Effect of provision on Pr(choose agreement)') + 
  scale_y_percent(limits = c(-.10, .10))
ggsave(here("repFile", "paper", "figures", "auth-inter-num.pdf"), 
       device = cairo_pdf)
ggsave(here("repFile", "JOP submission", "auth-inter-num.pdf"), 
       device = cairo_ps)
ggsave(here("repFile", "JOP submission", "auth-inter-num.pdf"), 
       device = cairo_pdf)


# sparse-reg plot results
load(here("repFile", "data", "sparse-reg-longrun.rda"))


# labels for plotting
labs2 = c('REMOVE', '
          Large land transfer', 
         'Small land transfer', 
         'Participate + seats', 
         '2 times min. salary', 
         'All rebels face prison', 
         'Rebels convicted of HR violations', 
         'Rebels pay victims', 
         'Rebel public apology', 
         'Air fumigation', 
         'REMOVE',
         'Child-rear X\n Participation + seats', 
         'Child-rear X\n 2 times min. salary', 
         'Child-rear X\n All rebels face prison', 
         'REMOVE')

# get output
pDat = 
  summary(s1)$table %>% 
  round(3) %>% 
  data.frame() %>% 
  mutate(label = rownames(.)) %>% 
  filter(`Posterior.Median` != 0) %>% 
  mutate(group = ifelse(stringr::str_detect(label, 'x'), 
                        'Interaction Effects', 
                        'Main Effects (AMCEs)')) %>% 
  mutate(lab2 = labs2) %>% 
  filter(lab2 != 'REMOVE') %>% 
  mutate(attribute = c('Reform', 
                       'Reform', 
                       'Elections', 
                       'Reintegration',
                       'Retributive justice', 
                       'Retributive justice', 
                       'Reparations', 
                       'Reparations', 
                       'Reform', 
                       'Elections', 
                       'Reintegration', 
                       'Retributive justice')) %>% 
  mutate(attribute = factor(attribute)) %>% 
  arrange(attribute) %>% 
  mutate(order = 1:nrow(.))


# throw out trailing zero
no_zero <- function(x) {
  y <- sprintf('%.1f',x)
  y[x > 0 & x < 1] <- sprintf('.%s',x[x > 0 & x < 1]*10)
  y[x == 0] <- '0'
  y[x > -1 & x < 0] <- sprintf('-.%s',x[x > -1 & x < 0]*-10)
  y
}

ggplot(pDat) + 
  geom_pointrange(aes(x = reorder(lab2,order), 
                      y = `Posterior.Median`,
                      ymin = `X5.`, ymax = `X95.`, 
                      colour = attribute, 
                      #shape = attribute, 
                      NULL 
                      )) + 
  coord_flip() + 
  facet_grid(rows = vars(group), scales = 'free') + 
  theme_light(base_family = "Roboto Condensed") + 
  geom_hline(yintercept =  0, lty = 2) + 
  labs(x = '', y = '', 
       # title = 'Results: LASSOPlus', 
       # subtitle = 'Estimated effects for non-zero coefficients.', 
       NULL) + 
  scale_color_manual(values = wesanderson::wes_palette(n = 5, 
                                                       name = 'Darjeeling1')) + 
  theme(legend.position = 'bottom',
        legend.title = element_blank(), 
        panel.border = element_blank(),
        plot.title = element_text(face = 'bold', size = 18),
        title = element_text(hjust = 0),
        axis.ticks = element_blank(),
        strip.background = element_rect(fill = 'Black'),
        strip.text = element_text(size = 12),
        axis.text.y = element_text(face = 'bold'))
ggsave(here("repFile", "paper", "figures", "lassoplus-long.pdf"), 
       device = cairo_pdf)
ggsave(here("repFile", "JOP submission", "lassoplus-long.eps"), 
       device = cairo_ps)
ggsave(here("repFile", "JOP submission", "lassoplus-long.pdf"), 
       device = cairo_pdf)
